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ABSTRACT 

Finite  element  modeling  is  being  adopted  in  the  design  of  ultrasonic  transducers  and  imaging  arrays. 
Impetus  is  accelerated  product  design  cycles  and  the  need  to  push  the  technology.  Existing  designs  are 
being  optimized  and  new  concepts  are  being  explored.  This  recent  acceptance  follows  the  convergence  of 
improvements  on  many  fronts:  necessary  computer  resources  are  more  accessible,  lean,  specialized 
algorithms  replacing  general-purpose  approaches,  and  better  material  characterization 

The  basics  of  the  finite  element  method  (FEM)  for  the  coupled  piezoelectric-acoustic  problem  are 
reviewed.  We  contrast  different  FEM  formulations  and  discuss  the  implications  of  each:  time-domain 
versus  frequency  domain,  implicit  versus  explicit  algorithms,  linear  versus  nonlinear.  Beyond  discussions 
of  the  theoretical  underpinnings  of  numerical  methods,  the  paper  also  examines  other  modeling  ingredients 
such  as  discretization,  material  attenuation,  boundary  conditions,  farfield  extrapolation,  and  electric 
circuits. 

Particular  emphasis  is  placed  on  material  characterization,  and  this  is  discussed  through  an  actual  "model- 
build-test"  validation  sequence,  undertaken  recently.  Some  applications  are  also  discussed. 
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1.  INTRODUCTION 

Until  recently,  medical  transducer  designers  relied  almost  exclusively  on  ID  analytical  models  and 
experimental  prototypes.  Now,  many  employ  comprehensive  finite  element  simulations  for  transient,  2D 
and  3D  analyses.  Effectiveness  of  their  modeling  is  proportional  to  accuracy  and  completeness  of  material 
measurements,  fidelity  of  geometrical  and  manufacturing  process  details,  and  the  modeler's  skill  with 
numerical  experiments  and  design  strategies.  Modeling  skills  rely  on  an  intuitive  understanding  of  the 
transducer's  operational  parameters  and  the  overall  design  problem.  This  requires  either  years  of 
experience  or  focused,  practical  training.  Without  such  skills,  modeling  is  often  used  poorly  and  results 
may  be  misinterpreted  or  erroneous,  leading  to  wasted  resources  and  market  opportunities.  As  personal 
computers,  modern  finite  element  algorithms,  and  user  interfaces  make  2D/3D  modeling  more  accessible, 
the  required  skill  set  needs  to  be  taught  quickly  and  efficiently  to  an  eclectic  mix  of  users. 

Therefore,  in  our  role  as  code  developers,  users,  and  trainers,  we  have  sought  a  simple  yet  practical  basis 
for  teaching  transducer  fundamentals  and  the  finite  element  modeling  paradigm.  The  most  direct  basis 
would  be  closed-form  solutions,  but  transducers  generally  have  too  many  characteristic  lengths  for  simple 
mathematical  analysis.  However,  note  that  each  segment  of  a  typical  medical  transducer  has  a  fairly  short 
extensional  resonance  length,  e.g.,  «  A/2  for  the  piezoceramic  and  ~  A/4  for  the  matching  layer(s).  These 
fractional-wave  dimensions  suggest  that  relatively  low  frequency  coupled  oscillators  rather  than 
propagating  waves  can  provide  a  simple,  intuitive  model  of  device  physics.  Although  oscillator  models 
cannot  be  quantitative,  in  general,  they  provide  the  simplest,  complete  representation  of  ID 
electromechanical  transducer  behavior.  This  approach  appears  to  have  been  passed  over  in  the  transducer 
literature  and  as  a  teaching  aid.  To  this  end.  Section  2  examines  simple  spring-mass  models  that 
demonstrate  basic  device  behavior  and  as  well  as  introduce  finite  element  fundamentals,  given  that  springs 
and  masses  are  the  simplest  "elements"  of  discrete  numerical  modeling. 
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We  then  proceed  in  Section  3  to  generalization  of  these  electromechanical  models,  leading  to  the  full  3D 
finite  element  formalism.  This  was  first  developed  by  Allik  and  Hughes  in  1970  [1],  Although  in  use  since 
then  for  analysis  of  low-frequency  underwater  projectors  [2],  its  adoption  in  the  medical  ultrasound 
community  remained  limited  until  the  early  1990's.  The  investigation  of  the  then  new  1-3  piezocomposites 
accentuated  the  need  for  such  comprehensive  modeling,  as  exemplified  by  the  work  of  Hossack  and 
Hayward  [3],  The  work  of  Lerch  [4]  emphasized  the  need  for  transient  response  modeling  and  non-uniform 
damping  in  realistic  applications,  features  lacking  in  then  available  commercial  software.  It  is  the  adoption 
of  explicit  wave  propagation  algorithms  in  PZFlex  [5-11]  however  that  made  realistic  transducer 
simulations  practical.  With  demonstrated  speed/size  advantage  factors  of  100  over  conventional  implicit 
algorithms,  broadband  imaging  transducer  models  became  tractable  on  desktop  computers.  The  point  is  that 
software  used  in  a  production  capacity  must  rely  on  specialized  algorithms  rather  than  general-purpose 
ones.  To  this  end,  we  contrast  different  FEM  formulations  and  discuss  the  implications  of  each:  time- 
domain  versus  frequency  domain,  implicit  versus  explicit  algorithms,  and  linear  versus  nonlinear  schemes. 

Beyond  discussions  of  the  theoretical  underpinnings  of  numerical  methods,  we  examine  in  Section  4 
modeling  issues  from  the  analyst's  perspective:  discretization,  material  dissipation,  boundary  conditions, 
farfield  extrapolation,  and  electric  circuits.  The  "quality"  and  "cost  effectiveness"  of  the  model  depends  on 
the  proper  use  of  these  ingredients.  We  review  their  underlying  assumptions  and  provide  guidelines  for 
optimal  use. 

Accuracy  of  the  finite  element  model  also  hinges  on  accuracy  of  the  material  constitutive  properties,  and 
those  provided  in  manufacturer  specification  sheets  are  often  incomplete,  if  not  inaccurate.  Section  5 
examines  material  characterization  issues  in  the  context  of  a  recently  undertaken  validation  exercise  based 
on  an  incremental  "model-build-test"  analysis  of  a  nonproprietary  ID  biomedical  imaging  array.  This 
sequence  of  component  and  device  validations  provides  an  excellent  opportunity  to  identify 
characterization  procedures  that  work  and  those  that  require  further  refinement,  all  while  displaying  its 
impact  on  the  resultant  finite  element  solution.  Ultimately,  validated  and  established  characterization 
protocols  are  critical  if  numerical  modeling  is  to  serve  as  a  "virtual  prototyping"  tool. 

Some  recent  applications  and  studies  are  briefly  discussed  in  Section  6,  mainly  to  highlight  some  of  the 
important  points  made  earlier,  to  mention  new  ones  that  could  not  be  discussed  in  detail  in  the  body  of  the 
paper,  and  to  provide  some  references  to  a  broader  range  of  applications. 


2.  COUPLED  OSCILLATOR  MODELS  OL  RESONANT  TRANSDUCERS 

Each  segment  of  a  typical  medical  transducer  has  a  fairly  short  extensional  resonance  length,  e.g.,  ~  X/2  for 
the  piezoceramic  and  «  X/4  for  the  matching  layer(s).  These  fractional-wave  dimensions  suggest  that 
relatively  low  frequency  coupled  oscillators  rather  than  propagating  waves  can  provide  a  simple,  intuitive 
model  of  device  physics.  Although  oscillator  models  cannot  be  quantitative,  in  general,  they  provide  the 
simplest,  complete  representation  of  ID  electromechanical  transducer  behavior.  As  such,  they  are  valuable 
in  illustrating  the  "fundamentals". 

To  this  end  we  examine  what  amounts  to  spring-mass  models  of  transducer  stack  elements.  These  are  the 
most  fundamental  electromechanical  analogs  and  can  be  developed  intuitively,  as  follows,  or  from  more 
rigorous  matrix  structural  analysis  and  finite  element  concepts.  In  particular,  piezoelectric  constitutive 
relations  for  the  spring  are  derived  from  the  continuum  relations  and  applied  to  the  simple  harmonic 
oscillator.  This  is  generalized  to  a  coupled  oscillator  representing  the  lowest-order  resonances  in  a 
piezoelectric  transducer  stack,  i.e.,  backing,  piezoelectric,  matching  layer(s),  and  water  load. 

2.1.  Equations  of  motion  and  constitutive  relations 

The  basis  for  approximate  dynamic  models,  in  general,  is  (1)  assumption  of  a  specific  strain  function,  e.g., 
constant  or  linear,  and  (2)  mass  distribution,  i.e.,  consistent  or  lumped.  This  is  the  foundation  of  matrix 
structural  analysis  [12],  the  predecessor  of  finite  element  methods.  The  strain  function  assumption  yields 
ordinary  differential  equations  in  time  by  reducing  the  infinite  degrees  of  freedom  to  a  finite  number,  while 
the  mass  distribution  assumption  simplifies  time  integration. 
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Consider  a  ID  piezoelectric  bar  with  length  L  and  area  A,  electroded  on  the  ends  and  poled  lengthwise. 
Longitudinal  displacement  is  u(x,t),  governed  by  the  partial  differential  equation 


.  c~u  cT 
pA — —  =  — 

A2  Ax 


(2.1) 


where  x  is  the  space  coordinate,  t  is  time,  p  is  material  density,  and  T(x,t)  is  longitudinal  stress.  Boundary 
conditions  are  specified  on  u  or  T,  along  with  initial  conditions.  The  piezoelectric  constitutive  relations 
between  stress  T,  strain  S  =  An/ Ax,  electric  field  E,  and  electric  displacement  D  are 

T=cES-eE  ,  D  =  ssE  +  eS  (2.2) 

where  cE  is  elastic  stiffness  under  constant  electric  field,  e  is  the  piezoelectric  stress  constant,  and  ss  is 
electric  permittivity  under  constant  strain,  e.g.,  see  [13].  Because  V-D  =  0  (divergence  condition  on 
electric  displacement)  and  the  absence  of  free  charge,  D(t)  is  uniform  over  the  bar. 


The  simplest  approximation  of  the  bar's  dynamics  follows  from  the  constant  strain  assumption  and  mass 
lumping.  This  is  equivalent  to  a  linear  spring  separating  equal  end  masses  m  =  pAL/2,  as  shown  in  Fig.  1 . 
End  displacements  uj  and  Lp  are  the  degrees  of  freedom,  whence  governing  equation  (2.1)  simplifies  to  the 
ordinary  differential  equations 


d2u,  d2u1 

m — —  =  F  +Fl  ,  m- 


dt~ 


dt~ 


=  ~Fspr  +  F2 


(2.3) 


where  Fspr  is  the  spring  force  and  Fj,  F2  are  external  forces  applied  at  each  end.  Multiplying  (2.2)  by  A 
and  defining  spring  force  Fsp,=-AT ,  electrode  charge  Q  =  AD,  spring  compression  u  =  u/  -  i<2  =  LS,  and 
voltage  V  =  LE  yields  the  piezoelectric  spring  constitutive  relations 

Fspr  ~  ~kEu~  dmu  +  eV  ,  Q=CSV  +eu  (2.4) 

kE=AcEIL  ,  e  =  Ae  I L  ,  Cs=Aes/L  (2.5) 

where  kE ,  e ,  and  Cs  are  spring  stiffness,  piezoelectric  force  constant,  and  capacitance,  respectively.  Note 
that  the  spring  equation  is  generalized  above  to  include  rate  dependent  mechanical  damping  through  dm. 
The  sign  convention  on  u  is  positive  for  compression  and  negative  for  extension.  The  spring  idealization, 
governing  equations,  and  constitutive  relations  are  illustrated  in  Figure  1 . 


T  =  cES  -  eE  ,  D  =  £ SE  +  eS  F  =  -kE(uj- u?)  +  eV  ,  Q  =  CSV  +  e(uj-U2) 


A  d2u  8T 
PA — r  =  — 

dt 2  dx 


U1  2—  u2 


Fig.  1.  Longitudinal  piezoelectric  oscillator  (left)  of  length  L  and  area  A,  and  the 
equivalent  simple  harmonic  oscillator  model  (right),  showing  the  corresponding 
piezoelectric  constitutive  relations  (above)  and  governing  equations  (below). 


2.2.  Simple  harmonic  oscillator  solutions 

Combining  the  spring-mass  equations  of  motion  and  the  spring  constitutive  relations  yields  two  ordinary 
differential  equations  governing  the  piezoelectric  oscillator.  For  example,  consider  the  symmetric  problem, 
i.e.,  F 1  -  -F2  =  FeXf  and  uj  =  -i<2,  whence  (2.3)  and  (2.4)  yield 

miil+2dmu1 +2kEul  =  eV  +  Fexl  ,  Q  =  CSV +2eu{  (2.6) 

Two  canonical  cases  of  interest  are  the  "send"  problem  defined  by  voltage  input  and  displacement  output, 
and  the  "receive"  problem  defined  by  external  force  input  and  voltage  output. 
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For  the  send  problem,  the  differential  equation  and  constitutive  relation  for  current  i  =  dQ /dt  are 

mul+2d,nul+2kEul=eV  ,  i  =  CsV+2eul  (2.7) 

When  the  prescribed  voltage  is  time  harmonic,  V{t)  =V0elwl ,  steady  solutions  are  ul{t)  =  uQei°”  and 
substituting  gives 


^o_  = _ £ _ 

V0  -mco2  +  j2dmco  +  2kE 

Electromechanical  response  is  given  by  the  electrical  impedance,  i.e.,  voltage  divided  by  current, 

(  r\i  E  __2  \ 


V_ 

i 


J 


1 


co  C  +2 eu0/V0 


J 


E  2 

2k  -  mco 


co  Cb(2kE  -ma)  +  2e 


for  d,„  =  0 


(2.8) 


(2.9) 


For  the  receive  problem  assume  that  current  vanishes,  hence,  Q  =  Qo  =  0  (without  loss  of  generality,  since 
Qo  produces  static  compression  and  voltage  that  are  removable  by  redefining  equilibrium  position  and 
ground  potential).  The  governing  equation  and  constitutive  relation  for  voltage  become 

miil  +2  dmul  +2k°ul  =  Fext  ,  CSV  +eu{  =0 
where  kD  is  the  stiffened  spring  constant 

kD=kE(l  +  K2)  ,  K2  = 

kECs  cEss 

and  K  is  the  piezoelectric  coupling  constant.  For  time  harmonic  external  force,  Fexl(t)  =  F0ejmt ,  steady 


(2.10) 

(2.11) 


solutions  are  u1(t)  =  u0e]a  ,  V(t)=VQeJ  and  substituting  gives 


1 


1 


F0  2k  -  mco  +  j2d  mco  2  k  -mco 

Output  sensitivity  is  defined  as  voltage  divided  by  force, 
{  ^ 


for  d,„  =  0 


V 


2e  m0 


-2e 


{  Cs (2k°  -  mco2) 


for  d,„  =  0 


J 


(2.14) 


(2.15) 


These  simple,  closed-form  solutions  exhibit  most  of  the  piezoelectric  resonator  characteristics  that  are  of 
interest  to  the  designer  and  modeler.  Of  course  they  are  crude  approximations  despite  being  complete 
representations  of  ID  transducer  physics.  In  contrast,  more  comprehensive  models  of  the  type  described  in 
this  paper  provide  more  complete,  quantitative  answers  but  without  simple  functional  relations. 


2.3.  Coupled  harmonic  oscillator  solutions 

Generalizing  the  single  piezoelectric  oscillator  described  above  to  a  coupled  oscillator  representing  a 
transducer  stack  is  straightforward  and  illustrated  in  Figure  2.  The  piezoelectric  ceramic  and  matching 
layer(s)  are  replaced  by  piezoelectric  or  elastic  springs  with  half  of  the  mass  lumped  at  each  end.  The 
water  and  backing  loads  are  represented  by  dashpots  with  coefficient  d  =  ImA  where  Im  is  mechanical 
impedance  of  the  load  medium  and  A  is  cross-sectional  area  of  the  stack.  This  follows  from  ID  wave 
theory,  which  states  that  pressure  p  and  velocity  x  are  related  as  p  =  -Im  x  where  Im  =  pc  is  the 
mechanical  impedance  of  the  medium. 

Consider  the  single  matching  layer  case.  Coefficients  and  spring  forces  are 
ml=\LlplAl  ,  m2=m1+m3  ,  m3=\L2p2A2 

db  =  PA’Lb\  .  dw  =  p„yLwA2  (2A4) 

Fl=—kf(ul—u2)  +  eyi  ,  F2  =  -k2(u2  - m3) 

where  vLb  and  vLw  are  wave  speeds  (longitudinal)  in  the  backing  and  water  load.  For  simplicity  we  ignore 
intrinsic  damping  in  the  matching  layer  and  piezoceramic.  Equations  of  motion  for  the  three  masses  are 
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(2.15) 


mliil  =Fl—  dbul  =  —kf(ul  -u2)  +  elVl  -dbitl 
m2u2  =  F2  —  Fl  =  —k2(u-,  -u3)  +  kf  (iq  —  «2)  —  e{V 
m3u3  =~F2  —  dn,u3  =  k2(u2  —u3)-dwii3 

The  input  is  driving  voltage  V  =  V0eim> ,  hence,  we  seek  the  time-harmonic  solutions 

rq  =  u^1  ,  u2  =  u02e]  ,  u3  =  Uq3€^ 

Substituting  yields  the  symmetric,  tri-diagonal  system  of  equations 


-  mlco 2  +  kf  +  jcodb 

-kf 

o  A 

fu  > 

u0l 

rn 

1 

—  m2co~  +  k2  +  kf 

—  k2 

u02 

II 

-i 

0 

V 

—  k2 

- m3co 2  +k2+  jcodw 

X<)3  ) 

,0, 

This  system  of  equations  representing  the  forced  vibration  problem  is  easily  solved  as 
uQ2  _  e^m^coi1  -  jcodb)(m3co~  —k2  —  jcodw ) 

V0  (»qnr  - k f  -  j cod b)(m2co2  —k2  — kf)(m3a> 1  —k2  —  jcodw)  - 

[kUm^o2  -  kf  -  j  cod  b)  -  (kf  )2  (m3co2  -  k2  -  j  cod  „,)  J 

e  +kE- T  -k  ^ 

uoi  _ _ qo _  ho3  _ _ qo _ 

Vo  —ni\Cd2  +  kf  +  jcodb  V0  -m3co2  +k2  +  jcodw 
Intrinsic  material  damping  in  the  piezoelectric  and  matching  layer  is  included  by  letting 
kf  — >  kf  +  jcodml  ,  k2  — »  k2  +  jcodm2 

in  these  equations,  where  dm\  and  dm2  are  the  mechanical  spring  damping  coefficients. 


(2.16) 


(2.17) 


(2.18) 


(2.19) 

(2.20) 


backing 


water 


sN|  db 


matching 

piezoceramic  layer(s) 


m2  h I  \U  \U 

m4) - 

/| — ►  1*2 

/) — -  Uj 

/j — ►  u4 

" — ^  k2 

d3 

V 

% 

H.' fn,3  J 

dr^u2 

\\ 

X 


Fig.  2.  Coupled  harmonic  oscillators  representing  the  basic  transducer  stack  with  a 
single  or  double  matching  layer. 
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On  this  basis,  closed-form  solutions  of  coupled  piezoelectric  oscillator  models  are  easily  studied.  They  are 
the  manifest,  low  order  limit  of  representations  ranging  from  electrical  analogs  like  Mason  and  KLM 
models,  to  ID  wave  propagation  models,  to  broadband  (transient)  2D  and  3D  finite  element  models.  In 
terms  of  utility,  such  models  have  proven  convenient  as  a  "trivial"  mathematical  basis  for  demonstrating 
generic  device  behavior,  and  more  specifically,  for  teaching  fundamentals,  illustrating  matching  layer 
design  issues,  and  qualitative  interpretations  of  transducer  experiments  and  numerical  simulations.  More  to 
the  point,  they  are  the  archetype  for  2D/3D  finite  element  transducer  models  and  an  intuitive  guide  to 
proper  electromechanical  device  modeling,  particularly  if  the  engineer's  background  has  focused  on 
electrical  analog  interpretations. 


3.  ALGORITHMIC  FORMULATION  OF  PIEZOELECTRIC  FINITE  ELEMENTS 

The  concepts  illustrated  by  the  simple  oscillator  models  above  are  generalized  to  full  3D  behavior  next, 
leading  to  the  complete  finite  element  formulation.  This  entails  a  dimensionality  extension  of  course,  but 
mainly  raises  the  issue  of  effective  computer-based  solution  schemes. 

3.1.  Governing  differential  equations  (strong  form) 

Conventional  polarized  ferroelectric  ceramics  used  in  the  manufacture  of  ultrasonic  transducers  are 
governed  by  the  constitutive  relations  for  linear  piezoelectricity  [13]  and  by  the  equations  of  mechanical 
and  electrical  balance.  The  electric  balance  is  assumed  instantaneous  and  decoupled,  i.e.,  the  quasi-static 
approximation  of  the  electric  field  relative  to  the  mechanical  field.  The  governing  equations  are  thus 
expressed: 


T  =  C£  S-er  E  ,  D  =  e  S  +  e'  E 

-  constitutive  equations- 

(3.1) 

pa  =  V  T 

-  momentum  balance  - 

(3.2) 

V  •  D  =  0 

-  electric  balance  - 

(3.3) 

with 

S  =  Vs  •  u  ,  E  =  -V<f> 

where  T,  S,  E,  D  are  the  mechanical  stress,  mechanical  strain,  electric  field  and  electric  displacement 
vectors,  respectively.  CE,  8s,  e  are  the  matrices  of  stiffness  constants  at  constant  electric  field,  of  dielectric 
constants  at  constant  strain,  and  of  piezoelectric  coupling  constants,  respectively,  u  is  the  mechanical 

displacement  vector  and  ii  =  <32u/<3?2  is  the  acceleration,  each  superposed  dot  "."  denoting  one  time 

differentiation.  (j>  is  the  electric  potential  (voltage).  To  complete  the  description  of  the  problem,  equations 
(3.1)-(3.3)  are  complemented  by  appropriate  boundary  conditions,  such  as  prescribed  displacements  or 
voltages,  and  applied  forces  or  electric  charge.  These  equations  are  the  3D  counterparts  to  (2.1)  and  (2.2), 
with  scalar  quantities  now  replaced  by  matrices. 

3.2.  Semidiscrete  finite  element  equations 

The  finite  element  method  (FEM)  is  an  approximation  of  the  governing  equations  that  is  particularly  well 
suited  to  computation.  Whereas  the  differential  form  of  the  governing  equations  requires  the  solution  to  be 
exact  at  every  point  in  space,  the  FEM  is  based  on  an  equivalent  variational  or  "weak"  statement  that 
enforces  the  "exactness"  of  the  solution  in  a  weighted  average  sense  over  small  sub-regions  of  the  space 
(the  finite  elements).  For  the  class  of  formulations  considered  here,  convergence  and  uniqueness  of  the 
solution  can  be  established  mathematically.  In  other  words,  error  bounds  on  the  approximation  can  always 
be  determined  and  the  approximation  can  always  be  improved  such  that  the  weighted  error  tends  to  zero  at 
the  limit,  or  equivalently,  the  finite  element  solution  tends  to  the  exact  solution.  In  practice  though,  one 
does  not  require  a  zero  error,  but  an  error  small  enough  to  be  insignificant  compared  to  other  sources  of 
uncertainty  (e.g.,  experimental  errors  in  determining  material  properties,  geometric  dimensions, 
manufacturing  tolerances),  commonly  lumped  under  the  "umbrella"  of  noise. 

The  FEM  requires  the  domain  of  the  problem  to  be  subdivided  into  small  discrete  finite  elements:  4-node 
quadrilateral  in  2D  or  8-node  hexahedron  in  3D,  for  example.  The  solution  sought  is  expressed  in 
polynomial  expansions  with  the  coefficients  of  the  polynomial  being  the  value  of  the  solution  field  at  the 
finite  element  nodes.  In  other  words,  the  FEM  solution  vector  consists  of  the  displacement  values  Uj  and 
electric  potential  values  (j);  at  nodes  i;  the  displacement  and  voltage  fields  at  arbitrary  locations  within 
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elements  are  determined  by  a  linear  combination  of  polynomial  interpolation  for  shape)  functions  Nu  and 
N^,  respectively,  and  the  nodal  values  of  these  fields  as  coefficients: 

u(x,  y,  z,t )  =  N„  ( x ,  y,  z).u(f)  =  N[  (x,  y,  z).ue  it)  (3.4a) 

y,z,t )  =  N^fx,  y,z).O(0  =  N^(x,  y,z).<Ds(f)  (3.4b) 

where  superscript  "e"  denotes  quantities  associated  with  a  given  element.  Equations  (3.4)  highlight  key 
characteristics  of  the  FEM  as  opposed  to  Rayleigh-Ritz  methods,  namely: 

1)  The  nodal  "unknowns"  of  the  problem  have  physical  significance  (e.g.,  displacement)  and  are  not  just 
expansion  coefficients. 

2)  FEM  interpolation  functions  are  local  or  element-based,  implying  the  solution  within  an  element  is 
entirely  determined  by  the  solution  at  that  element's  nodes.  It  is  this  localization  that  permits  element- 
by-element  operations,  and  therefore  allows  the  FEM  to  solve  large-scale  complex  problems  as  an 
assembly  of  tractable,  elemental  contributions. 

When  the  shape  functions  are  taken  to  be  linear,  the  strain  distribution  within  the  element  is  constant,  just 
like  the  spring  in  Section  2.  In  that  spirit,  the  quadrilateral/hexahedron  continuum  elements  can  be  viewed 
as  2D/3D  spring  and  mass  combinations. 


Incorporation  of  the  spatial  discretization  (3.4)  into  the  above  mentioned  variational  statement  results  in  a 
semidiscrete  finite  element  system  of  linear  algebraic  equations,  expressed  in  matrix  form  as  follows: 

M>  +  C>  +  K„u  +  =  F  (3.5a) 

K>  +  K#<X>  =  Q  (3.5b) 

where 


M 


uu 


A 


jpNfKdV 


-  mechanical  mass  matrix  -  (3.6a) 


K„„=Aj(VN:)rC£(VN:)£/yf 

=  A  J(VN;,)rer(VN;)c/ye 

e  Ve 

K#=Aj(VN;)V(VN;)jye 

Cuu  is  the  mechanical  damping  matrix,  F  and  Q  are  the  nodal  mechanical  force  and  electric  charge  vectors, 
respectively,  and  u  and  O  are  the  nodal  displacement  and  potential  vectors,  respectively.  The  scheme  by 
which  elemental  contributions  are  assembled  to  form  the  global  system  matrices  is  represented  by  the 

element  assembly  operator  A .  This  element  assembly  process  is  akin  to  the  simple  one  shown  in  (2.14), 
leading  to  (3.5)  which  is  the  3D  time-domain  counterpart  of  (2.17). 


-  mechanical  stiffness  matrix  -  (3.6b) 

-  piezoelectric  coupling  matrix  -  (3.6c) 

-  dielectric  stiffness  matrix  -  (3.6d) 


Equation  (3.5a)  governs  the  mechanical  or  elastic  portion  of  the  problem,  while  equation  (3.5b)  describes 
the  electrical  field,  and  both  are  coupled  through  the  piezoelectric  coupling  matrix.  For  passive  materials, 
the  coupling  is  null  and  equation  (3.5a)  fully  describes  the  behavior  of  elastic  materials.  We  note  that 
inviscid  and  irrotational  fluid  (acoustic)  media  are  sometimes  more  conveniently  described  by  potential- 
based  formulations  [4,5],  which  we  shall  not  describe  here  for  the  sake  of  brevity.  Equations  (3.5)  are 
referred  to  as  the  semidiscrete  FE  equations  in  that  space  has  been  discretized  whereas  time  is  still 
represented  as  a  continuous  function.  To  solve  such  a  system,  assumptions  and/or  approximations  on  the 
time  dimension  must  be  made. 


3.3.  The  time  dimension  and  associated  solution  algorithms 

Frequency-domain  analysis:  When  the  dynamic  phenomenon  is  steady-state,  with  periodic  forcing  function 
and  response  at  circular  frequency  0)  =  27if,  time  dependence  can  be  eliminated  from  the  problem  and  the 
system  unknowns  convert  to  harmonic  complex  variables: 

u  =  ueM,  ®  =  8Q/dt  =  >(.),  d2(.)/dt2  =  -®2(.)  (3.7) 
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The  FE  equations  (3.5)  then  reduce  to  a  complex  symmetric  non-hermitian  matrix  system  requiring  an 
implicit  solver.  Direct  implicit  solution  by  Gaussian  elimination  is  only  practical  in  2D  because  3D  leads  to 
prohibitively  large  system  bandwidth  and  memory  needs.  For  larger  problems,  iterative  solvers  are 
indicated.  For  problems  free  of  material  and  radiation  damping,  the  system  is  positive  definite  and  the 
conjugate  gradient  (CG)  method  is  appropriate.  In  realistic  situations  though,  various  attenuation 
mechanisms  (e.g.,  water  loading)  result  in  a  typically  indefinite  system  requiring  more  general  iterative 
solvers  such  as  GMRES  [14]  and  QMR  [15].  In  practice,  the  utility  of  frequency-domain  formulations 
diminishes  as  the  phenomenon  of  interest  involves  multi-modal  behavior,  especially  at  higher  frequencies. 
Even  in  transducers  intended  for  steady-state  operation,  one  typically  analyses  the  spectrum  around  the 
narrow  operational  band  to  insure  modal  "purity",  and  that  requires  resolving  the  response  over  several 
discrete  frequencies,  with  one  full  system  solution  for  each.  That  is  not  to  say  that  the  utility  of  examining 
the  data  in  the  frequency-domain  is  diminished,  since  it  is  a  concise  and  convenient  visualization  of 
complex  behavior,  but  rather  argues  to  the  inefficiency  of  the  algorithmic  approach  in  such  instances.  In 
others  words,  the  solution  domain  should  be  dictated  by  computational  efficiency  since  the  data  can  always 
be  viewed  in  any  desired  domain  through  relatively  speedy  post-processing  conversion  by  FFT.  Many  of 
the  early  FE  implementations  for  piezoelectric  dynamics  were  formulated  in  the  frequency  domain,  most 
likely  because  of  the  early  concentration  on  low-frequency  sonar  applications  and  the  relative  simplicity  of 
extending  real  arithmetic  elastostatic  FEM  to  complex  arithmetic  harmonic  elastodynamics  FEM. 

Eigenvalue/Eigenmode  extraction :  Classical  eigenanalysis  for  extraction  of  natural  frequencies  and 
associated  modal  shapes  also  requires  matrix  factorization,  whether  direct  or  iterative,  and  these  can  be 
found  in  many  textbooks.  It  should  be  noted  though  that  eigenanalysis  becomes  computationally  difficult  as 
modal  separation  diminishes  at  higher  resonances,  and  that  the  eigensolution  only  pertains  to  energy 
conserving  (i.e.,  undamped)  systems.  Attenuation  effects  due  to  water  loading  or  material  damping  are  not 
accounted  for.  When  damped  modes  are  of  interest,  one  needs  to  analyze  the  forced  vibration  problem  (in 
either  the  frequency  or  time  domain)  for  a  range  of  the  spectrum,  identify  resonances,  and  then  extract  the 
displacement  field  (i.e.,  mode  shape)  at  these  resonant  frequencies. 

Time-domain  analysis:  When  transient  or  broadband  signals  are  of  principal  interest,  the  temporal  evolution 
of  the  system  is  best  resolved  through  step-by-step  time  integration  schemes.  It  is  also  the  only  solution 
approach  if  nonlinear  phenomena  are  involved.  There  are  many  ways  to  determine  the  current  solution  at 
time  tn+i  from  known  solutions  at  the  previous  time  step  tn  (algorithms  involving  higher  order  time 
approximations  will  involves  several  past  time  levels,  but  these  tend  to  be  reserved  for  special  situations  in 
view  of  the  associated  computational  burden).  The  Newmark  family  of  time  integrators  is  widely  used  in 
mechanical  dynamics.  It  assumes  a  constant  acceleration  a  =  92u/9r  over  a  small  time  interval  (time-step) 
At  and  insures  2nd  order  accuracy  in  the  temporal  approximation.  It  takes  on  the  following  form  when 
applied  to  the  mechanical  FE  equation  (3.5a): 

W+i  =v„  +  Ma„  +a„+1)/2 

u„+i  =  u„  +  A  rv„  +  At 2  [(l  -  /?)a„  +  ]/2 

a„+i  =  _  M,™  +  —  c,™  +  K™  x  F„+i  +  C„„(^„  +  —  a„  j  +  K„„  u„  +  A rv„  +  -^-(l  -  Py„ 

v - v - J 

LHS 

F„  ,  =F„n  -K„,  0„„  (3.8) 

Different  choices  of  the  Newmark  parameter  /?  result  in  temporal  integrators  optimized  for  different  classes 
of  problems: 

Implicit  methods  (e.g.,  J3  =  1/4)  couple  current  solution  vectors,  hence,  the  global  system  of  equations  must 
be  solved  at  each  time  step.  The  LHS  in  (3.8)  involves  a  matrix  factorization  (e.g.,  by  Gaussian 
elimination)  that  is  an  expensive  operation  requiring  a  computational  effort  of  order  ()(n2„„,ie)  ~  0(n:' node) 
and  storage  of  order  O (n2node).  The  advantage  of  implicit  schemes  is  unconditional  stability  with  respect  to 
the  time  step.  Implicit  methods  are  typically  indicated  for  statics  (e.g.,  electrostatics  (3.5b)),  low-frequency 
mechanical  or  inertial  dynamics,  and  diffusive  processes  (e.g.,  thermal  diffusion)  described  by  parabolic 
PDEs,  where  temporal  gradients  are  substantially  smaller  than  spatial  ones. 
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Explicit  methods  decouple  the  current  solution  vectors  and  eliminate  the  global  system  solve,  but  they  are 
only  conditionally  stable,  i.e.,  there  is  a  time  step  limit  (CFL  condition  [16])  beyond  which  the  algorithm 
becomes  unstable.  In  elastodynamics,  the  CFL  time  step  limit  corresponds  to  the  shortest  transit  time  across 
any  element  in  the  mesh  (Afstab=min(h/c),  h  is  the  element  size  or  nodal  distance,  and  c  is  the  wavespeed). 
In  wave  phenomena,  the  desired  resolution  and  accuracy  require  a  time  step  smaller  than  one-tenth  the 
period  of  the  highest  frequency  of  interest,  a  requirement  no  less  stringent  than  that  imposed  by  the  CFL 
condition,  and  thus  removing  the  principal  advantage  of  implicit  methods.  It  is  this  limit  on  the  time  step, 
and  therefore  on  the  distance  traveled  during  each  interval,  that  allows  the  nodal  fields  to  decouple 
momentarily  during  a  time  step.  Although  not  usually  implemented  in  this  fashion  and  for  the  purpose  of 
this  discussion,  consider  an  explicit  scheme  as  a  Newmark  integrator  with  ft  =  0  and  Muu,  Cuu  diagonalized 
by  nodal  lumping.  In  this  case,  each  equation  in  system  (3.8)  can  be  integrated  independently,  i.e.,  in  a 
decoupled  fashion.  The  coupling  is  then  effectively  accounted  for  through  the  forces  on  the  right  hand  side 
of  the  equation,  and  these  are  vector  operations  that  are  much  less  demanding  in  computer  resources  than 
the  matrix  factorization  required  by  implicit  schemes.  In  actual  implementation,  the  matrices  shown  in  (3.8) 
are  not  assembled  and  stored,  since  an  explicit  scheme  naturally  structures  itself  into  a  series  of  element-by- 
element  operations  involving  global  vectors  only.  As  such,  storage  and  solution  requirements  scale  linearly 
with  the  number  of  nodes  in  explicit  schemes,  which  is  still  a  requirement  for  useful  3D  modeling. 
Element-by-element  operations  also  offer  a  major  advantage  in  code  parallelization,  which  is  gaining 
mainstream  interest  with  the  availability  of  multiple  processor  PCs. 

The  piezoelectric  problem  couples  a  mechanical  dynamics  process  (3.5a)  and  an  electrostatic  one  (3.5b). 
Computational  efficiency  would  suggest  that  a  mixed  explicit/implicit  scheme  is  optimal  for  the 
mechanical/electrical  problem.  This  mixed  scheme  has  been  demonstrated  to  achieve  a  2  orders  of 
magnitude  efficiency  gain  (computational  speedup  for  a  given  model,  or  model  size  for  a  given 
computational  time)  with  PZFlex  [5]  compared  to  conventional,  fully  implicit  FEM  implementations.  This 
advantage  is  expected  to  be  maintained  as  long  as  the  electromechanical  problem  maintains  a  wave 
propagation  characteristic,  and  increases  markedly  in  large-scale  applications.  To  illustrate  the  point, 
consider  the  naval  flextensional  transducer  depicted  in  Fig.  3.  This  high  power  PMN-driven  flextensional  is 
one  of  12  sources  constituting  a  towed  array,  and  the  corresponding  model  requires  about  10  million  finite 
elements.  An  analysis  of  such  large-scale  and  involving  material  nonlinearities  (PMN)  could  only  be 
undertaken  today  with  an  explicit  scheme.  Similar  large-scale  problems  exist  in  the  biomedical  arena  [10]. 


Fig.  3.  Behavior  of  a  single  flextensional  model  driven  in  water,  showing  radiating  wave 
pattern  (left).  View  of  mode  shape  in  the  single  flextensional  model  driven  at  4.5 
kHz  (right).  A  detailed  description  is  available  in  [9] 
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4.  FINITE  ELEMENT  MODELING  ISSUES  IN  TRANSDUCERS  AND  ARRAYS 

As  in  any  analysis  effort,  the  "quality"  of  the  answers  obtained  depends  not  only  on  a  good  choice  of 
methods  but  also  on  the  input  parameters  used  to  define  the  model.  Chief  among  these  is  material 
constitutive  properties,  which  we  discuss  in  detail  in  Section  5.  Other  parameters  control  the  level  of  model 
fidelity  to  the  physical  behavior,  typically  optimized  against  computational  cost.  It  is  often  argued  that 
finite  element  simulation  could  be  turned  into  a  "black  box"  if  computational  resources  were  infinite.  While 
this  statement  does  point  out  a  legitimate  constraint,  it  ignores  that  one  often  has  to  analyze  preliminary 
designs  with  tentative  data  and  that  different  levels  of  approximations  are  appropriate  for  different 
questions.  In  what  follows,  we  discuss  some  of  the  basic  modeling  issues  underlying  finite  element 
simulation  of  transducers  and  arrays,  and  provide  guidelines  for  the  best  use  of  these  representations. 

4.1.  Spatial  and  temporal  discretization 

Meshing  a  finite  element  model  defines  the  solution's  resolution.  For  wave  propagation  problems,  the 
discretization  must  resolve  the  shortest  wavelength  (i.e.,  highest  frequency)  of  interest.  This  is  analogous  to 
crystal  lattice  theory  [17]  where  the  periodic  discrete  structure  defines  a  low-pass  filter  with  the  highest 
frequency  propagated  (cutoff  frequency)  having  a  wavelength  equal  to  2h  (h=  finite  element  size  or  nodal 
distance);  this  is  also  known  as  the  Nyquist  limit  or  saw-tooth  pattern.  Obviously,  one  desires  adequate 
rather  than  borderline  resolution  of  the  frequencies  of  interest,  and  this  ranges  from  7,/h  =  8  to  20  finite 
elements  per  wavelength.  Discretization,  by  virtue  of  being  on  approximation  process,  introduces  a 
numerical  error  that  displays  a  non-physical  frequency  and  direction  dependent  dispersive  character.  A 
wave  pattern  is  well  resolved  with  7,/h  =  10  elements  per  wavelength  (Fig.  4)  and  this  produces  a  numerical 
dispersion  error  of  about  3%.  When  long  propagation  distances  are  involved  and  wavefront  or  pulse 
distortions  must  be  limited,  gridding  with  20  elements  per  wavelength  limits  the  numerical  error  to  less 
than  1%  [18-20].  Because  these  frequencies  propagate  throughout  the  model,  the  same  degree  of 
discretization  should  be  maintained  throughout  the  mesh  to  avoid  directionality  and  spurious  internal 
reflectors.  In  practice,  the  finite  element  aspect  ratio  should  remain  close  to  one,  but  can  safely  reach  2  or  3 
to  accommodate  geometric  constraints.  In  contrast,  a  diffusive  process  such  as  heating  tends  to  display 
localized  spatial  gradients  and  the  degree  of  mesh  refinement  may  vary  from  region  to  region  accordingly 
(Fig.  4). 


Fig.  4.  Spatial  discretization  for  wave  propagation  (left)  vs.  diffusion  (right)  problems. 


(a)  In  a  ID,  2nd  order  accurate  grid  (FE  or  FD). 


(b)  In  a  Pseudospectral  grid  integrated  in  time  by 
leapfrog,  Runge/Kutta,  and  Adams-Bashforth. 


Fig.  5.  Comparison  of  wavelet  time  histories  after  propagating  300  wavelengths  [10]. 


10 


Furthermore,  because  of  the  duality  of  space  and  time  in  the  wave  equation,  the  same  criterion  used  for 
spatial  discretization  applies  to  time  discretization,  and  means  that  the  optimal  time-step  should  be  at  the 
CFL  stability  limit.  In  practice  though,  only  90-95%  of  CFL  in  linear  cases  and  80%  of  CFL  in  nonlinear 
cases  are  achievable,  which  is  sufficient  for  transducer  analysis.  However,  when  long-range  wave 
propagation  must  be  accurately  described,  in  bioacoustic  models  for  example,  standard  linear  finite 
elements  are  no  longer  effective  because  of  the  accumulation  of  dispersion  errors.  Then,  schemes  such  as 
the  pseudospectral  method  using  higher-order  approximations  must  be  used  [10],  as  shown  in  Fig.  5. 

Frequency  content  is  not  always  the  controlling  factor  in  vibration  problems.  Geometric  features  and 
boundaries  sometimes  impose  more  stringent  requirements  on  discretization.  A  prime  example  of  this 
occurs  in  the  case  of  the  PZT  half  wavelength  thickness  resonator  in  a  medical  transducer  stack.  If 
discretization  is  set  by  wavelength  considerations  alone,  the  bar  width  would  be  divided  into  one  or  two 
finite  elements  across,  which  misses  the  lateral  stress  gradients  due  to  Poisson  effects.  That  lateral 
distribution  influences  coupling  with  an  adjacent  polymer  matrix,  for  example,  and  gridding  the  width  with 
6  nodes  is  suggested  if  such  effects  need  to  be  resolved. 

4.2.  Material  Damping 

Frequency-dependent  material  damping  is  an  important  issue  in  transducer  modeling  because  of  the  many 
polymers  involved  including  backing,  matching  layers,  composite  matrix,  and  lenses.  Damping  not  only 
affects  the  acoustic  signals  but  also  generates  heat.  In  the  frequency-domain,  the  damping  level  can  be 
specified  at  each  discrete  frequency,  without  concern  as  to  overall  frequency  power  laws  and  wide 
spectrum  characterization.  Classical  time-domain  finite  element  damping  models  are  chosen  for  their 
operational  and  algorithmic  properties  rather  than  their  phenomenological  behavior.  These  models  are 
restricted  in  their  frequency  dependence,  but  can  vary  from  one  element  to  another: 

[a]  Mass  proportional  damping  yields  an  attenuation  per  unit  distance  that  is  constant  with  frequency, 
a  =  aff*1).  The  equivalent  critical  damping,  which  is  a  measure  of  attenuation  per  wavelength  expressed 
as  a  percent  of  critical  (level  required  for  full  attenuation  over  1  wavelength),  follows  an  inverse 
frequency  dependence  t  cc  1/f, 

[b]  Stiffness-proportional  damping  yields  a  quadratic  dependence  of  attenuation  a  =  a(f2)  on  frequency 
(critical  damping  E,  <x  f), 

[c]  Rayleigh  damping  allows  for  a  linear  combination  of  the  mass  and  stiffness  proportional  models. 

[d]  Three -parameter  viscoelastic  models  provide  a  slightly  more  complex  behavior,  with  the  attenuation 
exponent  varying  from  2  to  0  with  increasing  frequency,  but  are  still  amenable  to  FEM  treatment. 


These  models  provide  an  adequate  fit,  in  general,  over  a  range  of  frequencies  as  shown  in  Fig.  6.  The 
Rayleigh  damping  model  is  versatile  in  that  it  provides  a  2  parameter  fit.  Substantial  improvement  in 
matching  experimental  results  is  obtained  when  damping  properties  are  specified  independently  for  the 
volumetric  (or  longitudinal)  and  deviatoric  (or  shear)  components,  for  most  materials.  Polymers  and 
rubbers  for  example  exhibit  much  greater  attenuation  of  their  shear  components. 


s.o 

6.0 

4.0 

2.0 

0.0 


Classical  viscous  absorption  models  Rayleigh  damping  fits  to  f 1  ■ 1  and  f 1  5  power  laws 


Fig.  6.  Frequency  dependence  of  various  damping/attenuation  models  used  in  time 
domain  wave  propagation  calculations  (left),  and  Rayleigh  damping  fits  to 
frequency  power  laws.  All  models  are  made  to  match  a  prescribed  value  at  5 
MHz. 
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Viscoelastic  models  with  more  general  frequency  power  laws  can  be  formulated  [21],  but  carry  a  high 
computational  overhead  because  of  the  retarded  integrals  involved.  Research  into  alternative  forms 
amenable  to  explicit  time -domain  calculations  [22]  is  ongoing  and  warrants  further  attention,  but  urgency 
on  that  front  is  tempered  by  the  current  limitations  of  experimental  characterization  of  polymers. 

4.3.  Boundary  conditions 

It  is  often  impractical  to  model  the  full  extent  of  the  transduction  device  and  the  surrounding  acoustic 
media.  In  fact,  considering  that  typical  ultrasonic  wavelengths  are  in  the  millimeter  range  at  best,  high 
resolution  finite  element  simulations  would  be  computationally  exorbitant  unless  the  domain  was  restricted 
to  the  area  of  interest.  Truncation  by  an  artificial  boundary  is  required  for  domains  large  compared  to  the 
characteristic  wavelength.  Appropriate  boundary  conditions  need  to  be  imposed  on  the  truncation  boundary 
to  simulate  the  behavior  of  a  "continuing"  medium.  Waves  incidents  on  the  boundary  need  to  exit  the 
computational  domain  with  no  spurious  reflections,  consistent  with  the  fact  that  the  boundary  is  a 
mathematical  construct  and  not  an  impedance  mismatch  zone  (Fig.  7).  Such  boundary  conditions  have  been 
termed  radiation,  transmitting,  absorbing,  non-reflecting  or  silent  conditions.  Reviews,  extended  reference 
lists,  or  comparative  studies  on  this  subject  can  be  found  in  [23-25], 


Absorbing  BC 


Backing  Layer 

Fig.  7.  Finite  computational  domain  bounded  by  an  artificial  truncation  boundary 
(dashed  line)  at  which  an  absorbing  boundary  condition  is  applied. 

There  are  no  exact  absorbing  conditions  applicable  to  all  situations:  frequency-domain,  time-domain, 
nonlinear  wave  propagation.  In  the  linear  case,  retarded  potential  integral  formulations  define  an  exact 
boundary  condition,  but  one  which  is  spatially  and  temporally  non-local.  Non-local  schemes  couple  all 
degrees-of-freedom  at  the  truncation  boundary  and  all  past  time  steps,  which  is  impractical  in  realistic  finite 
element  studies  because  of  their  expense.  Computationally  attractive  local  schemes  are  based  on  varying 
degrees  of  approximation,  with  assumptions  on  the  spatial  decay  rate  of  the  radiating  wave,  its  angle  of 
incidence  on  the  truncation  boundary,  and  its  wave  speed.  The  common  feature  to  all  local  radiation 
conditions  is  that  they  are  asymptotically  exact  at  high  frequency,  i.e.,  the  wavelength  is  shorter  than  the 
scale  of  the  boundary.  The  simplest  and  least  accurate  of  these  is  the  water  load  impedance  condition,  p  =  - 
I m  x  where  Im  =  pc  used  in  (2.14),  which  is  exact  in  ID  but  only  asymptotically  exact  in  the  high- 
frequency/farfield  limit  in  2D/3D.  The  only  currently  available  boundary  treatment  that  is  equally 
applicable  to  linear  and  nonlinear  wave  propagation  is  one  recently  proposed  by  Sandler  [26]  and 
implemented  in  PZFlex,  for  which  we  coin  the  acronym  MINT  (Material  Independent  Non-reflecting 
Treatment)  condition.  A  brief  derivation  of  the  MINT  condition  is  given  below. 


The  normally  propagating  part  on  a  wave  traveling  towards  the  truncation  boundary  is  given  by  the 
Hadamard  identity  for  outgoing  waves  dldn=  (-1  /cn)d/dt,  where  cn  is  the  unknown  wave  phase  velocity  in 
the  direction  of  the  outward  normal  n  to  the  boundary.  The  Hadamard  identity  applies  to  the  traction  vector 
x  =  T.n  at  the  boundary,  A„x/A n  =  -  (l/c„)(A,x/Af),  so  that  the  equation  of  motion  or  momentum  balance 
yields  a  change  in  boundary  nodal  velocity: 

=  =  (4.i) 

At  p  An  pcn  At 
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(4.2) 


The  Hadamard  identity  also  applies  to  the  boundary  velocity  field  such  that: 

_  1  V 

An  cn  At 

Combining  (4.1)  and  (4.2)  by  eliminating  the  unknown  wavespeed  c„  yields  the  MINT  condition: 


This  absorbing  boundary  condition  makes  no  assumptions  regarding  material  constitutive  properties  (i.e., 
nonlinearities  can  be  accommodated),  the  boundary  geometry,  or  angle  of  incidence  of  the  scattered  wave. 
The  MINT  condition  has  been  shown  to  perform  as  well  as  a  4th  order  paraxial  absorber  [27],  with  lower 
computational  overhead  and  less  impact  on  stability. 

Because  of  its  unusual  accuracy  and  its  natural  fit  within  a  finite  element  approach,  we  should  note  the 
radiation  boundary  condition  recently  developed  by  Berenger,  and  coined  the  Perfectly  Matched  Layer  or 
PML  [28].  Although  originally  developed  for  electromagnetics,  it  has  been  shown  to  be  applicable  and 
particularly  effective  in  acoustic  wave  propagation  [10,29].  The  PML  construct  is  not  universally 
applicable  (e.g.,  nonlinear  acoustics  cannot  be  accommodated  and  the  elastodynamic  case  has  not  been 
formulated  to  date),  but  it  fills  a  crucial  niche  in  the  feasibility  of  large-scale  bioacoustic  models. 

In  practice,  one  needs  to  exercise  some  care  in  the  placement  of  absorbing  boundary  conditions.  They 
should  be  placed  at  some  distance  away  from  active  sources  where  wave  patterns  are  complex  and 
gradients  severe.  This  also  avoids  near  grazing  incidence  on  the  boundary,  where  accuracy  generally 
degrades.  When  cross-talk  and  coupling  effects  are  of  interest,  the  wave  path  must  be  kept  in  the 
computational  domain,  since  a  wave  absorbed  by  the  non-reflecting  boundary  condition  cannot  be 
reintroduced  at  another  point.  Finally,  and  when  in  doubt,  a  useful  sanity  check  consists  of  plotting  the 
pressure  or  stress  waves  patterns  in  the  entire  mesh  and  verifying  that  spurious  reflections  from  the 
boundary  are  at  second  order  error  levels  and  not  physically  meaningful  levels.  In  such  circumstances, 
animations  serve  not  only  to  explain  complex  interaction  phenomena  but  also  to  validate  the  model. 

Absorbing  treatments  are  not  the  only  boundary  conditions  that  afford  model  reduction.  Symmetry  and 
periodicity  concepts  are  equally  crucial  in  effective  solution  approaches.  Approximating  finite  spatial 
periodicity  by  an  infinite  periodicity  is  often  used  with  the  caveat  that  edge  effects  are  filtered  out. 
Similarly,  reduction  of  a  3D  structure  with  a  "long"  dimension  to  a  2D  plane-strain  model  is  also  a  routine 
modeling  approach.  Even  when  the  aspect  ratio  does  not  justify  such  a  reduction,  a  2D  analysis  still  serves 
the  purpose  of  a  "first  cut"  that  can  guide  more  exhaustive  and  expensive  3D  studies.  Examples  in  Section  6 
and  the  listed  references  make  evident  the  extensive  and  effective  use  of  such  reductions. 


4.4.  Farfield  and  nearfield  extrapolation 

Results  at  some  distance  away  from  the  transducer  are  often  of  interest,  such  as  beam  patterns,  at  pulse- 
echo  reflectors  or  focal  points.  Solution  by  finite  elements  would  require  field  calculations  in  the 
intervening  region  between  the  source  and  the  distant  output  points,  which  is  of  impractical.  A  better 
alternative  is  offered  by  exterior  integral  formulations  that  only  require  the  discretization  of  a  surface 
bounding  the  "source"  region  and  only  calculate  the  solution  at  specified  output  points.  As  such,  they  can 
be  viewed  as  extrapolation  methods.  Integral  formulations  exist  for  propagation  through  homogeneous 
elastic  media  and  even  for  multi-layered  elastic  media.  In  practice  though,  the  simplest  and  probably  most 
useful  integral  equations  describe  the  radiation  through  homogeneous  acoustic  media  (dilatational  waves 
only),  and  we  limit  the  discussion  here  to  these  instances. 


Time-domain  Kirchhoff  integral  equation:  the  exact  expression  for  the  time -dependent  pressure  at  any  point 
Xq  in  the  exterior  infinite  fluid  surrounding  the  vibrating  body  (transducer),  enclosed  by  a  surface  .S’,  is  given 
by  the  Kirchhoff  integral  equation: 


p(x0>t)=  Jj 


dv„(xs,r)  dG  ,  x 

p^r^G*^p(x-T) 


dr 


with  the  Green's  function  G  given  by: 


(4.4a) 
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-  in  2D  (4.4b) 

2 - 1  +  R/c )2  -  (fl/c)2 

where  R  denotes  the  distance  between  the  surface  sampling  points  xs  and  xo,  c  is  the  fluid  medium 
wavespeed.  The  pressure  in  the  field  is  obtained  from  this  convolution  integral  of  pressure  and  normal 
acceleration  dvjdt  time  histories  at  the  "sampling"  surface.  Since  the  Kirchhoff  equation  makes  no 
approximations  beyond  the  homogeneity  of  the  surrounding  unbounded  acoustic  medium,  the  solution  is 
valid  whether  the  field  point  is  in  the  nearfield  or  farfield. 

Although  the  integral  theorem  requires  S  to  be  a  closed  surface,  often  in  practice  we  take  it  to  be  a  plane 
surface  on  the  front  face  of  the  model  ignoring  contributions  from  the  other  sides.  This  modeling 
approximation  is  reasonable  to  the  extent  that  an  effective  aperture  can  be  defined,  but  caution  needs  to  be 
exercised  when  the  model  only  represents  a  portion  of  the  actual  radiating  device.  The  more  substantial  the 
sideway  energy  leakage  (e.g.,  shear  waves  in  matching  layers)  resulting  in  a  larger  aperture,  the  less  valid  is 
such  a  truncation.  Another  useful  simplification  in  analysing  medical  arrays  is  to  model  the  lens  and 
surrounding  water  as  one  homogenous  acoustic  medium.  Tthe  surface  data  is  then  sampled  at  an 
intermediary  level  between  the  top  of  the  matching  layer  and  the  absorbing  boundary. 

Frequency-domain  farfield  beam  patterns :  Beam  patterns  are  by  definition  regarded  as  frequency-domain 
angular  pressure  distributions  on  a  circle  in  the  farfield.  Although  the  Kirchhoff  equation,  or  its  frequency- 
domain  counterpart  known  as  the  Helmholtz  integral  equation,  could  be  used  to  calculate  beam  profiles, 
computational  expense  can  be  reduced  by  making  use  of  the  fact  that  /?—»«.  The  governing  integral 
equation  reduces  then  to  the  Rayleigh-Sommerfeld  diffraction  equation 

p.,  (ffl,  0)  =  J^e^cos (6>)f  e  'kRp(xs  )dS  in  2D  (4.5) 

which  can  be  further  simplified  with  assumptions  of  ray  acoustics  and  calculated  efficiently  using  a  spatial 
FFT  [30].  Expressions  similar  to  (4.5)  also  exist  for  the  3D  case.  Because  of  assumptions  inherent  in  this 
formulation,  the  sampling  surface  is  taken  to  be  a  plane,  and  concerns  similar  to  those  in  the  Kirchhoff 
case,  regarding  the  definition  and  extent  of  the  sampling  surface,  must  be  taken  into  account.  Finally,  it 
should  be  noted  that  the  definition  of  a  beam  pattern,  based  on  unimodal  single  frequency  assumptions, 
does  not  always  coincide  with  the  beam  experimentally  obtained  by  sensing  pulse  (as  opposed  to  CW) 
maxima  along  a  circle  in  the  farfield.  They  do  coincide  usually  when  the  transducer  response  is  dominated 
by  a  single  well  separated  resonance. 

4.5.  Electric  circuits 

Ultrasonic  transducers  are  invariably  connected  to  some  supporting  electronics,  frequently  through  a 
coaxial  cable.  During  the  design  process,  it  is  typically  necessary  to  model  the  drive  circuitry,  the  reception 
electronics  and  the  intervening  cable.  When  designing  the  electronics,  it  is  most  efficient  to  compute  the 
impulse  response  of  the  transducer  and  use  this  in  a  circuit  design  code.  When  looking  at  the  transducer 
details,  it  is  more  efficient  to  model  the  electronics  directly.  Fortunately,  a  few  lumped  parameter  circuit 
elements  (resistors,  inductors,  capacitors,  and  ideal  transformers)  usually  suffice.  Figure  8  shows  generic 
drive  (send)  and  receive  circuits. 

In  the  time  domain  approach,  electrical  boundary  conditions  on  the  electrodes  are  replaced  by  a  coupled  set 
of  equations  relating  the  voltage  and  charge  (or  their  time  derivatives)  on  the  electrodes  to  the  voltage  and 
charge  (or  their  time  derivatives)  throughout  the  circuit.  The  electrical  boundary  conditions  (open,  ground, 
applied  voltage  or  current)  then  apply  to  the  circuit  rather  than  to  the  electrode.  We  note  that  the  charges 
and  potentials  at  each  of  the  circuit  elements  are  fully  coupled  to  the  nodal  values  throughout  the  FE  model 
at  each  timestep.  Careful  algorithmic  design  is  required  to  solve  this  system  accurately  and  efficiently. 

Electrodes  themselves  are  typically  modeled  as  voltage  constraints,  i.e.,  the  voltage  on  each  electrode  node 
is  constrained  to  be  equipotential.  This  is  almost  always  an  accurate  approximation,  though  windowing 
could  be  accounted  for  if  need  be  and  so  can  unusual  electrical  resistance.  The  mechanical  effects  of 
electrodes  are  typically  neglected  too.  The  electrode  mass  and  stiffness  could  be  accounted  for,  but  they  are 
usually  negligible  when  compared  to  those  of  the  ceramic. 


a_S(r -t  +  R/c).m  3p. 
4  kR 
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(a)  Generic  drive  circuit  similar  to  Panametrics  pulser. 


(c)  Generic  receive  circuit. 


rON  cb  COAX 


Ron  =  "ON"  resistance  of  MOS  FET 
Rd  =  Damping  resistor  -  changes 
bandwidth  of  output  pulse 
Rin  =  Input  impedence  of  pre-amp 
Cg  =  Blocking  capacitor 
Cq  =  Static  capacitance  of  transducer 
T  =  Turn  on  tine  of  switching 
MOS  FET 


(b)  Equivalent  circuit  representation. 


Fig.  8.  Typical  (a)  send  circuit,  (b)  its  equivalent  series  circuit  representation  appropriate 
for  incorporation  in  a  finite  element  model,  and  (c)  typical  receive  circuit. 


5.  MATERIAL  CHARACTERIZATION  AND  INCREMENTAL  VALIDATION 

Before  large-scale  finite  element  models  can  be  effectively  utilized  for  the  design  of  ultrasonic  imaging 
systems,  it  is  essential  that  these  modeling  tools  be  fully  validated  and  key  implementation  issues 
identified.  Furthermore,  there  is  a  need  for  a  well-documented,  nonproprietary  "primer",  that  transducer 
designers  can  use  to  establish  confidence  in,  and  expertise  with,  the  time-domain  finite -element-modeling 
paradigm.  The  complete  primer  will  consider  issues  such  as  material  characterization,  model  configuration 
and  model  implementation  (e.g.  introduction  of  bondlines  and  electrodes). 

The  following  sections  describe  a  recent  PZFlex  validation  exercise  based  around  the  incremental 
“ model-build-tesf  ’  analysis  of  a  nonproprietary,  5MHz,  ID  biomedical  imaging  array  (Fig.  9).  This  type  of 
incremental  analysis  permits  transducer  models  to  be  developed  that  describe  the  transducer  at  each  distinct 
stage  of  the  manufacturing  process.  Consequently,  deviations  between  experimental  results  and  finite 
element  predictions  can  be  more  readily  analyzed  and  the  transducer  models  refined  as  appropriate,  i.e., 
sources  of  accumulative  error  are  minimized. 

5.1.  Material  Characterization 

Accuracy  of  finite-element  analysis  is  ultimately  dependent  upon  the  accuracy  of  the  dielectric, 
piezoelectric  and  elastic  properties  used  to  represent  the  model’s  constituent  materials.  Consequently,  in  the 
course  of  this  work,  particular  emphasis  is  placed  on  material  characterization  issues. 
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Fig.  9.  Diagram  of  validation  array  assembly  (2  PZT  slivers  per  physical  array  element) 

Piezoceramic  properties  are  obtained  via  the  application  of  curve -fitting  techniques  to  IEEE  standard 
resonator  measurements  [31].  The  material  properties  are  then  cross-checked  by  comparing  the 
experimental  impedance  responses  for  the  IEEE  standard  resonators  with  the  corresponding  PZFlex 
predictions.  In  the  work  described  here,  matching  layer  and  backing  properties  were  obtained  via  a 
combination  of  through-transmission  water  tank  measurements  [32]  and  measurements  made  with  a  pair  of 
contact  shear-wave  probes.  Unfortunately,  these  measurement  methods  have  several  distinct  drawbacks  and 
hence,  more  comprehensive  and  accurate  measurement  schemes  are  currently  being  developed  [33]. 

Piezoceramic  Characterization :  The  IEEE  standard  on  piezoelectricity  [34]  identifies  certain  geometrical 
shapes  that  may  be  used  to  facilitate  the  measurement  of  a  material’s  elastic,  dielectric  and  piezoelectric 
properties.  For  piezoelectric  materials,  there  are  5  standard  resonator  geometries  (Fig.  10).  These  resonator 
samples  are  specifically  designed  so  as  to  isolate  certain  types  of  resonant  behavior.  Consequently,  it  is 
possible  to  measure  those  material  properties  that  are  strongly  coupled  to  a  particular  resonant  mode.  The 
equations  used  to  determine  the  material  properties  (as  given  in  the  IEEE  piezoelectric  standard  [34])  have 
idealized  derivations,  and  assume  that  the  material  is  lossless.  In  practice,  all  real  materials  possess  certain 
loss  mechanisms,  and  hence  the  calculated  properties  will  be  subject  to  certain  inaccuracies.  A  refinement 
to  this  method  has  been  developed  by  researchers  at  the  Royal  Military  College  of  Canada  and  employs 
curve-fitting  techniques  to  more  accurately  determine  a  material’s  properties.  The  software  package,  PRAP 
[31],  was  used  to  perform  this  analysis  and  the  extracted  material  properties  for  Motorola  3203HD  PLZT 
are  given  in  reference  [35]. 


The  IEEE  resonators  shown  in  Figure  10  were  modeled  in  PZFlex  using  the  measured  properties  given  in 
[35],  Figures  1  la-1  Id  show  the  correlation  between  the  experimental  electrical  impedance  response  and 
the  corresponding  PZFlex  predictions  for  a  selection  of  these  resonator  samples.  In  all  cases,  PZFlex  is  seen 


(d)  Length  thickness  mode  (LTE)  (e)  Thickness  shear  mode  (TS) 

Fig.  10.  IEEE  standard  piezoelectric  resonators:(a)  TE,  (b)  LE,  (c)  RAD,  (d)  LTE,  (e)  TS. 
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(c)  Frequency  (MHz) 
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(d)  Frequency  (MHz) 


Fig.  11.  Impedance  magnitude  response  for  resonators:  (a)  TE  (upper  left),  (b)  LE  (upper 
right),  (c)  RAD  (lower  left),  and  (d)  LTE  (lower  right). 


to  demonstrate  excellent  correlation  with  experimental  results.  The  simulated  in-air  impedance  response 
for  the  thickness  extensional  mode  resonator  (TE)  correctly  predicts  the  spurious  modal  activity  lying 
between  the  electrical  resonance  at  4.1MHz  and  the  mechanical  resonance  at  4.7MHz  (Fig.  11a).  These 
spurious  resonances  are  due  to  lateral  modes  that  are  supported  by  the  resonator’s  physical  geometry.  If 
these  parasitic  modes  are  too  strongly  coupled  to  the  particular  resonance  of  interest,  the  extracted  material 
properties  will  be  inaccurate.  The  TE  resonator  used  in  the  current  work  effort  was  specifically  chosen  such 
that  its  dimensions  fully  satisfy  the  guidelines  laid  down  in  the  IEEE  piezoelectric  standard.  Since  the  other 
resonators  have  greater  modal  separation,  their  impedance  response  curves  appear  cleaner  and  more 
unimodal. 


It  is  important  to  note  that  a  given  set  of  IEEE  standard  resonators  will  provide  material  properties  that 
were  measured  over  a  range  of  different  frequencies.  For  example,  the  TE  resonator  provides  cjj  at 
4.7MHz  whereas  the  LTE  resonator  yields  d13  at  200kHz.  Typically,  all  material  properties  will  exhibit 
some  degree  of  frequency  dependence,  however,  with  care,  it  is  possible  to  select  a  set  of  properties  that 
give  consistent  results  over  a  wide  range  of  frequencies. 

Matching  Layer  and  Backing  Block  Characterization:  The  validation  array  considered  in  this  paper  has  a 
double  matching  layer  and  light  acoustic  backing  attached  to  its  upper  and  lower  surfaces  respectively  (see 
Figure  9).  Accurate  characterization  of  these  passive  materials  is  just  as  important  as  for  the  piezoceramic. 
The  material  properties  that  need  to  be  measured  are  longitudinal-wave  velocity  &  attenuation,  and  shear- 
wave  velocity  &  attenuation.  Longitudinal  velocity  and  attenuation  measurements  are  relatively 
straightforward  and  may  be  accomplished  via  a  simple  through-transmission  experiment.  Unfortunately, 
measurement  of  the  shear  properties  is  typically  much  more  difficult  and  potentially  less  accurate.  Shear 
properties  are  normally  obtained  by  attaching  a  pair  of  shear  wave  transducers  to  opposite  sides  of  a  test 
specimen  and  propagating  a  shear  wave  through  the  sample.  Unfortunately,  it  often  proves  difficult  to  get 
good  coupling  of  energy  between  the  transducers  and  the  sample.  Consequently,  the  measured  values  are 
subject  to  considerable  inaccuracies.  Furthermore,  these  measurements  are  typically  narrowband  so  results 
are  only  valid  over  a  narrow  range  of  frequencies.  Wu  at  the  University  of  Vermont  is  currently  refining  a 
wide-band,  through-transmission,  water  tank  characterization  technique  [33]  for  passive  isotropic  materials. 
It  provides  both  velocity  and  attenuation  data  over  a  wide  range  of  frequencies.  Until  material  properties 
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obtained  via  this  method  become  available,  the  characterization  techniques  described  in  [32]  will  be  used  to 
characterize  the  array’s  matching  layer  and  backing  block  materials  (see  reference  [35]  for  measured 
material  properties). 

5.2.  PZFlex  Validation  -  Incremental  Array  Samples 

Once  all  the  active  and  passive  materials  have  been  accurately  characterized,  it  is  feasible  to  proceed  with 
PZFlex  analysis  of  the  transducer  array  assembly  shown  in  Figure  9.  This  array  is  constructed  from 
Motorola  3203HD  PLZT  and  has  48  individual  elements,  each  comprising  two  sub-elements  as  shown  in 
Figure  9.  A  double  layered,  sub-diced  matching  layer  was  adopted  and  a  “%/zf”  acoustic  backing 
(Z«2.5Mrayl)  was  bonded  to  the  bottom  of  the  device  to  help  improve  bandwidth  characteristics  while 
maintaining  reasonable  transmit  sensitivity.  For  validation  purposes,  rather  than  attempting  to  model  the 
entire  array  assembly,  the  analysis  is  currently  restricted  to  the  set  of  incremental  array  components  shown 
in  Figure  12(a).  Each  sample  corresponds  to  half  an  individual  array  element,  and  was  fabricated  using  the 
same  techniques  used  to  construct  the  complete  array.  Consequently,  any  process  dependent  effects  should 
also  be  observed  in  the  incremental  samples. 

The  first  step  of  the  validation  exercise  is  to  confirm  that  the  PZFlex  prediction  for  each  of  these  sub-units 

agrees  with  experimental  values.  Figure  12(b)  shows  both  the  experimental  ( - )  and  simulated  (0 . 0) 

electrical  impedance  response  curves  for  Sample- 1.  At  various  stages  throughout  the  fabrication  process  the 
piezoceramic  slivers  are  subject  to  elevated  temperatures  and  other  process  effects.  These  conditions  cause 
the  piezoceramic  to  partially  depole,  i.e.  its  piezoelectric  coupling  constants  will  be  reduced.  Consequently, 


Fig.  12.  (a)  Diagram  of  incremental  array  components  with  sample  numbers  shown,  and 

electrical  impedance  response  for  (b)  Sample-1,  (c)  Sample-2,  and  (d)  Sample-3. 
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5%  depoling  was  assumed.  From  the  cross-plotted  results  shown  in  Figure  12(b),  excellent  correlation 
between  experiment  and  simulation  are  seen.  Figure  12(c)  shows  the  response  for  Sample-2  (single  320pm 

sliver  of  3203HD  with  a  160pm  inner  matching  layer).  In  this  figure,  the  curve  denoted  ( - )  shows  the 

experimental  result,  however,  there  are  2  simulated  responses.  The  first  curve,  ( - ),  was  obtained  using 

the  nominal  matching  layer  properties  given  in  reference  [35],  whereas  the  second  simulated  curve  (0 . 0) 

used  velocities  in  the  matching  layer  which  were  reduced  by  10%.  Furthermore,  a  10pm  bondline  was  also 
included  in  the  calculation.  Modified  properties  give  much  better  correlation  with  experiment  than  the 
original  nominal  values.  It  was  mentioned  earlier  that  the  matching  layer  properties  are  often  difficult  to 
measure  accurately  and  that  more  accurate  measurement  methods  are  currently  being  developed.  We  expect 
to  have  more  accurate  material  properties  in  the  near  future  and  anticipate  that  they  will  further  improve  the 
correlation  with  experiment.  Sample-3  adds  a  70pm  outer  matching  layer  to  the  Sample-2  configuration. 
The  impedance  characteristics  for  this  device  are  shown  in  Figure  12(d)  and  the  simulated  result  with 
bondlines  and  modified  material  properties  is  seen  to  compare  well  with  the  experimental  result.  It  is 
interesting  to  observe  that  small  changes  in  matching  layer  properties  (longitudinal  or  shear)  and  the 
thickness  of  the  bondline  can  have  considerable  impact  on  overall  device  response.  This  once  again 
emphasizes  the  requirement  for  rigorous  material  characterization  and  accurate  experimental 
measurements.  Samples  4  &  5  correspond  to  the  addition  of  the  backing  block  to  Samples  2  &  3 
respectively  (for  impedance  results  refer  to  reference  [35]). 


6.  APPLICATIONS:  VISUALIZATION  AND  VIRTUAL  PROTOTYPING 

A  picture  is  worth  a  thousand  words,  and  a  movie  or  animation  of  deformed  shapes,  pressure  fields,  electric 
fields,  etc.  is  often  worth  much  more.  Practically  speaking,  model  visualization  is  the  ability  to  “see”  the 
operation  of  a  “virtual”  device  at  convenient  scales  and  speeds.  It  provides  one  of  the  most  compelling 
reasons  for  numerical  modeling.  In  addition  to  its  technical  merits,  visualization,  provide  the  best  means  of 
communicating  concepts  and  designs  to  other  engineers,  customers,  and  management.  This  section 
describes  some  scenarios  and  applications  illustrating  the  process. 

6.1.  Composites 

Typically,  preliminary  transducer  design  is  based  on  experience  with  similar  devices,  rules-of-thumb,  and 
simplified  ID  models.  Once  a  preliminary  design  has  been  decided  on,  it  is  appropriate  to  build  and  run  a 
finite  element  model.  The  computed  impedance  versus  frequency  curve  will  often  reveal  some 
unanticipated  resonances  that  may  affect  device  performance  adversely.  In  this  case,  the  next  logical  step 
is  to  compute  and  display  or  animate  the  deformation  shapes  at  these  resonances.  Based  on  shape 
information  and  parameter  studies,  the  designer  is  usually  able  to  modify  the  design  in  order  to  minimize  or 
eliminate  the  spurious  modes. 

A  1-3  composite  transducer  example  from  a  matching  layer  study  [8]  is  presented  in  Figures  13-15.  The 
transducer  is  from  an  undersea  imaging  array  designed  by  Ultrex  Corp,  built  by  Materials  Systems  Inc.  and 
described  in  [36].  Fig.  13  shows  measured  and  calculated  impedance  curves  for  the  case  of  a  water-loaded 
diced  matching  layer  with  filled  and  unfilled  kerfs,  while  Fig.  14  shows  deformation  shapes  at  peaks  on  the 
beam  pressure  versus  frequency  plot  for  the  filled  case.  The  design  frequency  is  350  kHz.  An  unexpected 
mode  is  seen  in  Fig.  13  at  500  kHz  for  the  case  of  air  kerfs  (left  side).  The  extremes  of  deformation  at  this 
frequency  are  pictured  in  Fig.  15,  showing  a  highly  localized,  matching  layer  mode.  This  is  a  lateral  mode 
that  couples  strongly  to  bending  modes  of  the  PZT  pillars  in  the  1-3  composite.  Computer  studies  show 
that  this  mode  disappears  for  a  water  load  but  if  there  is  any  shear  stiffness  in  the  load  medium,  as  assumed 
in  the  calculation,  e.g.,  an  RTV  lens,  then  this  type  of  spurious  mode  can  be  supported. 

Several  finite  element  studies  of  composite  transducers  have  appeared  in  the  literature  in  the  recent  past, 
and  good  examples  of  1-3  connectivity  can  be  found  in  [3,37,38],  and  2-2  connectivity  in  [39], 
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Fig.  13.  Measured  and  simulated  electrical  impedance  for  the  water-loaded  MSI  coupon, 
with  and  without  matching  layer  kerf  filler. 


ninn  mil 

B:  360  kHz 


Fig.  14.  Response  charaterization  of  the 
composite  model,  showing  detailed 
spectrum  and  mode  shapes  for  the  1.6 
mm  thickness  matching  layer. 


Fig.  15.  Localized  matching/pillar  mode 
shapes  at  500  kHz  in  Fig.  13  for  air 
kerfs;  calculated  by  a  3D  model. 


6.2.  Therapeutics 

Therapeutic  ultrasound  is  another  application  where  finite  element  modeling  and  visualization  can  be  used 
to  great  advantage,  both  in  transducer  design  and  studies  of  ultrasound  fields  within  the  body.  Here,  the 
goal  is  to  deliver  a  high  intensity  ultrasonic  field  to  the  target  region  while  minimizing  its  effects  on 
surrounding  tissue.  As  an  example.  Figure  16  from  [7]  compares  pressure  fields  from  uniformly  driven  and 
phased  transducer  designs.  In  the  latter  case,  phasing  was  used  to  move  the  focal  region.  Results  were 
poor  due  to  strong  cross-talk  between  the  continuous  transducer  rings. 
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from  a  radially  poled,  2  mm  thick,  PZT- 
5H  spherical  cap  with  10  cm  aperture 
and  radius  (f/1).  Ten  equal-area  annular 
electrodes  are  driven  by  1.0  volt  at  1.0 
MHz,  uniformly  (left)  or  phased  (right). 


Fig.  17.  Axisymmetric  model  of  an  ocular 


tumor,  showing  focused  ultrasound  beam 
(above)  and  temperature  distribution  at  3 
sec  (below).  Note  focal  temperature  in 
excess  of  100  °C.  The  transducer  has  a  4 
cm  aperture  and  9  cm  focal  length. 


Because  of  high  intensities  at  the  focus,  material  nonlinearities  in  tissue  become  an  issue.  These  include 
both  compressive  nonlinearity  that  is  commonly  described  by  the  first  nonlinear  term  (B/A)  in  an  expansion 
of  the  pressure  density  relation,  and  cavitation  under  tensile  stresses.  Both  of  these  are  readily  modeled  in 
the  time  domain  [7],  We  note  that  material  nonlinearity  causes  the  generation  of  higher  harmonics. 
Damping  in  tissue  increases  with  increasing  frequency,  so  these  in  turn  induce  greater  losses  and 
consequently  higher  heat  generation.  At  extremely  high  intensities  such  as  those  found  in  lithotripsy,  or  at 
very  long  propagation  distances,  shock  phenomena  must  also  be  modeled. 


Transducer  and  wave  field  calculations  are  sometimes  the  means  rather  than  the  end  of  a  simulation  [7], 
Figure  17  shows  the  computed  pressure  field  (upper  half)  near  the  focus  of  a  therapeutic  ultrasound 
transducer.  From  this  we  calculate  the  thermal  heat  generation  per  unit  time  at  each  point  in  the  model  and 
then  solve  the  bioheat  equations  for  evolution  of  temperature  with  time.  Temperature  at  3  seconds  is  shown 
in  the  lower  half  of  the  figure.  High  temperatures  can  be  used  to  destroy  diseased  tissue,  but  collateral 
damage  to  healthy  tissue  by  thermal  or  cavitation  mechanisms  should  be  minimized.  Simulations  provide  a 
convenient  method  to  refine  treatment  strategies  prior  to  in  vivo  validations  on  laboratory  animals. 


6.3  Prototyping 

A  last  example  is  used  to  illustrate  virtual  prototyping  of  transducers.  Fig.  18  shows  a  3D  model  of  a 
Tonpilz  device  for  low  frequency  sensing  in  air.  This  classical  design  is  usually  used  for  water-loaded 
applications.  The  model  consists  of  a  tail-mass,  a  stack  of  four  PZT  rings,  and  a  conical  head-mass,  all  held 
together  by  a  tension  bolt  through  the  axis.  A  thick  matching  layer  with  stiffening  ring  is  bonded  to  the 
head  mass. 


One  issue  of  interest  for  nominally  axisymmetric  devices  is  bending  modes  caused  by  nonsymmetric 
influences.  One  common  influence  is  nonuniform  driving  by  the  piezoelectric  rings.  This  was  studied  by 
introducing  pseudorandom  variation  in  the  coupling  constant  around  each  ring.  Both  impedance  and 
velocity  spectra  were  examined  to  identify  nonsymmetric  modes.  The  principal  bending  mode  was  found 
at  7  kHz  and  shown  on  the  left  side  of  Fig.  18.  The  principal  longitudinal  mode  was  around  14  kHz,  shown 
on  the  right  side  of  Fig.  18.  Higher  harmonics  of  the  bending  mode  were  not  apparent.  Effects  of  mount 
locations  and  nonsymmetric  mounting  fixtures  were  also  studied.  Such  analyses  are  readily  done  with 
numerical  models  and  tend  to  reduce  the  need  for  an  exhaustive  set  of  prototype  devices.  Nonetheless, 
validation  against  experiment  is  mandatory. 
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Fig.  18.  Examples  of  bending  (left)  and  longitudinal  (right)  modes  of  an  air-coupled 
Tonpilz  transducer.  A  thick,  flexible  matching  layer  is  bonded  to  the  face  of  the 
conical  head-mass. 


7.  CONCLUSIONS 

This  paper  was  intended  as  a  primer  on  various  issues  pertinent  to  finite  element  modeling  of  ultrasonic 
transducers.  The  breadth  of  the  topic  is  vast  and  multi-disciplinary  in  nature,  and  this  overview  is  perforce 
selective.  We  have  attempted  to  review  the  basic  algorithmic  background  and  practical  modeling  issues  to 
the  extent  they  affect  simulation  strategies  and  capabilities.  We  have  also  tried  to  convey  the  idea  that 
modeling  is  not  an  exercise  independent  of  experimentation,  whether  to  determine  the  requisite  material 
properties  or  to  develop  intuitive  understanding  of  the  transducer's  operational  parameters.  Much  can  be 
done  to  improve  the  accessibility  of  advanced  simulation  techniques,  through  intuitive  graphical  interfaces, 
standardized  design  "templates",  and  robust  algorithms  that  require  less  interaction  between  the  end-user 
and  the  "numerics".  We  do  not  foresee  however  that  increased  accessibility  renders  obsolete  efforts  to 
develop  modeling  skills  and  basic  understanding  of  the  underlying  analytical  methods.  From  our  current 
perspective  at  least,  these  often  effectively  complement  design  skills:  a  good  model  reflects  a  good 
understanding  of  the  physics. 
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